FOR-LOOPS-ARCHIVE

5.1 Convolution

I have 2 types of data here, one being used as an example if they are actually equal or the function is working. While the other one is used for benchmarking.

# To elaborate an example
set.seed(120)
x <- rnorm(10, mean = 10, sd = 2)
y <- rexp(10, rate = 1.5)

# For Benchmarking
set.seed(120)
xx <- rnorm(5e3, mean = 10, sd = 2)
yy <- rexp(5e3, rate = 1.5)

#include <stdio.h>
#include <stdlib.h>

void convolve_c(double *x, double *y, int *n, int *m, double *result) {
    int result_size = *n + *m - 1;

    for (int i = 0; i < result_size; i++) {
        int j_min = (i >= *m - 1) ? i - (*m - 1) : 0;
        int j_max = (i < *n - 1) ? i : *n - 1;

        result[i] = 0;
        for (int j = j_min; j <= j_max; j++) {
            result[i] += x[j] * y[i - j];
        }
    }
}

convolve_c <- function(x, y) {
    totalSize <- length(x) + length(y) - 1
    result <- double(totalSize)  # Initialize result vector
    result <- .C("convolve_c", 
                 as.double(x), as.double(y), 
                 as.integer(length(x)), as.integer(length(y)),
                 result=double(totalSize))$result
    return(result)
}

convolve_c(x, y) |> head(10)

 [1]  1.866342 10.531917 12.357334 14.054101 34.520841 39.345087 41.662725
 [8] 55.926328 68.759844 65.767219

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector convolve_cpp(NumericVector x, NumericVector y) {
    int n = x.size();
    int m = y.size();
    int result_size = n + m - 1;
    NumericVector result(result_size);

    for (int i = 0; i < result_size; i++) {
        int j_min = (i >= m - 1) ? i - (m - 1) : 0;
        int j_max = (i < n - 1) ? i : n - 1;
        double sum = 0.0;
        // Unroll the inner loop for better performance
        for (int j = j_min; j <= j_max - 4; j += 4) {
            sum += x[j] * y[i - j] +
                   x[j + 1] * y[i - (j + 1)] +
                   x[j + 2] * y[i - (j + 2)] +
                   x[j + 3] * y[i - (j + 3)];
        }
        // Handle remaining elements if any
        for (int j = j_max - (j_max - j_min) % 4; j <= j_max; j++) {
            sum += x[j] * y[i - j];
        }
        result[i] = sum;
    }

    return result;
}

convolve_cpp(x, y) |> head(10)

 [1]  1.866342 10.531917 12.357334 14.054101 34.520841 39.345087 41.662725
 [8] 55.926328 68.759844 65.767219

function convolve_vector(x::Vector{Float64}, y::Vector{Float64})
    n = length(x)
    m = length(y)
    result = zeros(Float64, n + m - 1)
    paddedX = zeros(Float64, n + m - 1)
    paddedY = zeros(Float64, n + m - 1)

    paddedX[1:n] = x
    paddedY[1:m] = y

    for i in 1:(n + m - 1)
        for j in 1:i
            result[i] += paddedX[j] * paddedY[i - j + 1]
        end
    end

    return result
end

convolve_vector (generic function with 1 method)

convolve_jl <- JuliaCall::julia_eval("convolve_vector")
convolve_jl(x, y) |> head(10)

 [1]  1.866342 10.531917 12.357334 14.054101 34.520841 39.345087 41.662725
 [8] 55.926328 68.759844 65.767219

use extendr_api::prelude::*;

#[extendr]
fn convolve_rs(x: Vec<f64>, y: Vec<f64>) -> Vec<f64> {
    let n = x.len();
    let m = y.len();
    let result_size = n + m - 1;
    let mut result = vec![0.0; result_size];

    for i in 0..result_size {
        let j_min = if i >= m - 1 { i - (m - 1) } else { 0 };
        let j_max = if i < n - 1 { i } else { n - 1 };
        let mut sum = 0.0;
        
        let mut j = j_min;
        while j + 3 <= j_max {
            sum += x[j] * y[i - j]
                + x[j + 1] * y[i - (j + 1)]
                + x[j + 2] * y[i - (j + 2)]
                + x[j + 3] * y[i - (j + 3)];
            j += 4;
        }
        
        for j in j..=j_max {
            sum += x[j] * y[i - j];
        }
        
        result[i] = sum;
    }

    result
}

convolve_rs(x, y) |> head(10)

 [1]  1.866342 10.531917 12.357334 14.054101 34.520841 39.345087 41.662725
 [8] 55.926328 68.759844 65.767219

subroutine convolve_fortran(x, n, y, m, result)
    implicit none
    integer, intent(in) :: n, m
    real(8), intent(in) :: x(n), y(m)
    real(8), intent(out) :: result(n + m - 1)
    real(8) :: paddedX(n + m - 1), paddedY(n + m - 1)
    integer :: i, j

    ! Initialize paddedX and paddedY with zeros
    paddedX = 0.0d0
    paddedY = 0.0d0

    ! Copy elements of x and y to paddedX and paddedY, respectively
    paddedX(1:n) = x
    paddedY(1:m) = y

    ! Compute convolution
    do i = 1, n + m - 1
        do j = 1, i
            result(i) = result(i) + paddedX(j) * paddedY(i - j + 1)
        end do
    end do
end subroutine convolve_fortran

convolve_fortran <- function(x, y) {
  totalSize <- length(x) + length(y) - 1
  result <- .Fortran("convolve_fortran", 
                     as.double(x), as.integer(length(x)), 
                     as.double(y), as.integer(length(y)),
                     result=double(totalSize))$result
  return(result)
}

convolve_fortran(x, y) |> head(10)

 [1]  1.866342 10.531917 12.357334 14.054101 34.520841 39.345087 41.662725
 [8] 55.926328 68.759844 65.767219

convolve_r <- function(x, y) {
    n <- length(x)
    m <- length(y)
    result <- rep(0, n + m - 1)
    paddedX <- rep(0, n + m - 1)
    paddedY <- rep(0, n + m - 1)
    
    paddedX[1:n] <- x
    paddedY[1:m] <- y
    
    for (i in 1:(n + m - 1)) {
        for (j in 1:i) {
            result[i] <- result[i] + paddedX[j] * paddedY[i - j + 1]
        }
    }
    
    return(result)
}

convolve_r(x, y) |> head(10)

 [1]  1.866342 10.531917 12.357334 14.054101 34.520841 39.345087 41.662725
 [8] 55.926328 68.759844 65.767219

import numpy as np

def convolve_py(x, y):
    n = len(x)
    m = len(y)
    result = np.zeros(n + m - 1)
    paddedX = np.zeros(n + m - 1)
    paddedY = np.zeros(n + m - 1)
    
    paddedX[:n] = x
    paddedY[:m] = y
    
    for i in range(1, n + m):
        for j in range(1, i+1):
            result[i-1] += paddedX[j-1] * paddedY[i-j]
    
    return result

convolve_py <- reticulate::py$convolve_py
convolve_py(x, y) |> head(10)

 [1]  1.866342 10.531917 12.357334 14.054101 34.520841 39.345087 41.662725
 [8] 55.926328 68.759844 65.767219

Benchmark

bench::mark(
     C = convolve_c(xx, yy),
     `C++` = convolve_cpp(xx, yy),
     Julia = convolve_jl(xx, yy),
     Rust = convolve_rs(xx, yy),
     FORTRAN = convolve_fortran(xx, yy),
     R = convolve_r(xx, yy),
     Python = convolve_py(xx, yy),
     check = F
)

# A tibble: 7 × 6
  expression      min   median `itr/sec` mem_alloc `gc/sec`
  <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
1 C            40.1ms  41.48ms   23.8      356.2KB        0
2 C++         12.65ms  15.71ms   62.5       85.8KB        0
3 Julia        80.2ms  87.65ms   11.2       83.9KB        0
4 Rust       591.53ms 591.53ms    1.69        83KB        0
5 FORTRAN     75.31ms   75.8ms   13.2        274KB        0
6 R             9.21s    9.21s    0.109    273.6KB        0
7 Python       47.62s   47.62s    0.0210   206.2KB        0